**********************************************************************************************************
* The (Null) Effects of the Russian Invasion of Ukraine on Europeans' Attitudes Toward Democracy
* Enrique Hernández & Macarena Ares 
* Published at Research & Politics
* Last updated: September 2023
* Whole syntax must be "runned" at once 
* Original analyses conducted with Stata 16.1 MP and Macbook Pro 2.3 GHz Intel Core i5
**********************************************************************************************************

*+++++++++++++++++++++++++++++++PREAMBLE+++++++++++++++++++++++++++++++++++++* 

capture: cd "" // Set path to directory where replication package was unzipped

global original "00_original_data/"  
global do "01_do_files/"  
global output "03_output/"  
global aux "05_auxiliary_output/"    
global gtrends "06_gtrends/"

frames reset

set graphics off

******************************
* ADOs (uncomment to install)
******************************
*ssc install fre, replace
*ssc install estout, replace 
*ssc install blindschemes, replace 
*ssc install ebalance, replace
*ssc install coefplot, replace
*net install dm79.pkg,  from(http://www.stata.com/stb/stb56) 
*net install dm49.pkg,  from(http://www.stata.com/stb/stb39)
*ssc install mmat2tex, replace 
*ssc install rwolf2, replace
*ssc install addplot, replace

*****************************************
* Seeds and parameters for replications
*****************************************
local seed_corrections = 06061988
local replications = 250 // Set to 250 to replicate 

*****************************************************************
* Appending original ESS self-completiton and face-to-face data
*****************************************************************
use "${original}ESS10.dta", clear

append using "${original}ESS10SC.dta"

****************************************
* Merging with data from contact forms
****************************************
merge 1:1 idno cntry using "${original}ESS10CF.dta", gen(merge1)

merge 1:1 idno cntry using "${original}ESS10SCCF.dta", gen(merge2) update

drop if merge1 == 1 & merge2 == 1
drop if merge1 == 2
drop if merge2 == 2

***************************
* Formating date variables
***************************

* Face-to-face finish date
generate interview_face = dofc(inwde)
format %td interview_face

* Face-to-face start date
generate interview_face_start = dofc(inwds)
format %td interview_face

* Self-completion finish date
generate interview_selcomp = dofc(questcmp)
format %td interview_selcomp

* Self-completion start date
generate interview_selcomp_start = dofc(scwsds)
format %td interview_selcomp

*******************
* Naming countries
*******************
replace cntry = "Austria" if cntry == "AT"
replace cntry = "Bulgaria" if cntry == "BG"
replace cntry = "Switzerland" if cntry == "CH"
replace cntry = "Czech Republic" if cntry == "CZ"
replace cntry = "Germany" if cntry == "DE"
replace cntry = "Estonia" if cntry == "EE"
replace cntry = "Spain" if cntry == "ES"
replace cntry = "Finland" if cntry == "FI"
replace cntry = "France" if cntry == "FR"
replace cntry = "Greece" if cntry == "GR"
replace cntry = "Croatia" if cntry == "HR"
replace cntry = "Hungary" if cntry == "HU"
replace cntry = "Iceland" if cntry == "IS"
replace cntry = "Italy" if cntry == "IT"
replace cntry = "Lithuania" if cntry == "LT"
replace cntry = "Montenegro" if cntry == "ME"
replace cntry = "Macedonia" if cntry == "MK"
replace cntry = "Netherlands" if cntry == "NL"
replace cntry = "Norway" if cntry == "NO"
replace cntry = "Poland" if cntry == "PL"
replace cntry = "Portugal" if cntry == "PT"
replace cntry = "Serbia" if cntry == "RS"
replace cntry = "Sweden" if cntry == "SE"
replace cntry = "Slovenia" if cntry == "SI"
replace cntry = "Slovakia" if cntry == "SK"

******************
* Data management
******************

* Country numeric variable
encode cntry, gen(cntrynum)

* Activity variable for balancing weights
gen active = 0 
replace active = 1 if pdwrk == 1 | edctn == 1 | cmsrv == 1

* Data management (so that results fit nicely in plots)
gen age_10 = agea/10
gen eduyrs_10 = eduyrs / 10 

label var age_10 "Age (Years/10)"
label var eduyrs_10 "Education (Years/10)"
label var active "Active (working or studying)"

* Voted last election (for balance tests)
gen voted = vote 
recode voted (2/3=0)

tab vote voted, m

* Generating variable for mode of survey administration 
gen face = 1 if name == "ESS10e02_1" 
replace face = 0 if name == "ESS10SCe01"

* Generating variable attempts to survey (for balance tests). Only for face-to-face countries 
foreach var of varlist resulb*{
	gen rec_`var' = `var'
	recode rec_`var' (1=0) (2/8 = 1)
}

egen attempt_aux1 = rowtotal(rec_resulb1-rec_resulb37), missing

gen attempts = attempt_aux1 + 1

gen attempts_10 = attempts / 10 

* Generating variable refusals to participate. Oly face-to-face countries
foreach var of varlist outnic*{
	gen refu_`var' = `var'
	recode refu_`var' (1=0) (2=1) (3=0) (4=1) (5/13=0)
}

egen refusals = rowtotal(refu_outnic1-refu_outnic37), missing

replace refusals = 0 if refusals == . & face == 1

* Region numeric variable
drop if region == "999999"
encode region, gen(regionnum)

*****************************
* Shorter variable labels
*****************************
label var implvdm "Importance democracy"
label var accalaw "Strong leader above law"
label var fairelc "Free elections"
label var dfprtal " Partisan offer"
label var gptpelc "Retrospect. accountability"
label var cttresa "Rule of law"
label var rghmgpr "Minority rights"
label var medcrgv "Satisfaction democracy"
label var gvctzpv "Protection poverty"
label var grdfinc "Red. income differences"
label var votedir "Referenda"
label var viepol "Views people prevail over elite"
label var wpestop "Will of people cannot be stopped"
label var stfdem "Satisfaction with democracy"
label var keydec "Key decisions national not EU"
label var attempts_10 "Attempts to survey (N/10)"
label var refusals "Refusals to participate"
label var voted "Voted last national election"

**********************************
* Checking overlap with invasion
* Invasion date 24 February 2022 
**********************************

* Locals for date and initial bandwidth arround invasion date (1 month)
local invasion_date = date("24feb2022","DMY")
local invasion_ll = `invasion_date' - 30
local invasion_ul = `invasion_date' + 30 

* Face-to-face countries 
tab cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=.

* Self-completion countries
tab cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != . 

* Capturing countries with overlap for next analyses 
levelsof cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=. & cntry != "Iceland", local(cntries_face)

levelsof cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != ., local(cntries_self) 

******************************************************************
* Droping respondents who start in control and finish in treatment
******************************************************************

* Face to face countries
tab cntry if interview_face != interview_face_start & interview_face_start != .

* Self-administered countries (web based surveys only)
tab cntry if interview_selcomp != interview_selcomp_start & interview_selcomp_start != .

* Face-to-face countries
drop if interview_face_start  < `invasion_date' &  interview_face > `invasion_date' & interview_face_start != .

* Self-administered countries
drop if interview_selcomp_start < `invasion_date' &  interview_selcomp > `invasion_date' & interview_selcomp != .

***********************************************************
* Distribution of respondents arround invasion by country
***********************************************************

* Face-to-face countries

* Generating plots by country with the full sample
foreach x of local cntries_face{
	
	sum interview_face if interview_face >= `invasion_ll' & interview_face < `invasion_date' & interview_face !=. & cntry == "`x'"  // Cases before (30 days bandwidth)
	local before_bw = `r(N)'

	sum interview_face if interview_face <= `invasion_ul' & interview_face > `invasion_date' & interview_face !=. & cntry == "`x'"  // Cases after (30 days bandwidth)
	local after_bw = `r(N)'
	
	hist interview_face if cntry == "`x'"	,	///
											xaxis(1 2)	///
											freq	///
											xline(`invasion_date', lpattern(dash))	///
											xlabel(, format(%tdDD/NN/YY) axis(1) labsize(vsmall)) 	///
											xlabel(`invasion_date', format(%tdDD/NN/YY) axis(2) notick labsize(vsmall))	///
											xtitle("", axis(2))	///
											xtitle("Interview date", axis(1))	///
											ytitle("")	///
											ti("`x'", size(small))
	
	graph save "${aux}interviews_`x'", replace
}

* Generating plot with all the countries
graph combine	"${aux}interviews_Switzerland"	///
				"${aux}interviews_Greece"	///
				"${aux}interviews_Italy"	///
				"${aux}interviews_Montenegro"	///
				"${aux}interviews_Macedonia"	///
				"${aux}interviews_Netherlands"	///
				"${aux}interviews_Norway"	///
				"${aux}interviews_Portugal"	
				
graph export "${output}distribution_interviews_ftfc.eps", replace


* Self-completion countries 

* Generating plots by country with the full sample
foreach x of local cntries_self{
	
	sum interview_selcomp if interview_selcomp >= `invasion_ll' & interview_selcomp < `invasion_date' & interview_selcomp !=. & cntry == "`x'"  // Cases before (30 days bandwidth)
	local before_bw = `r(N)'

	sum interview_selcomp if interview_selcomp <= `invasion_ul' & interview_selcomp > `invasion_date' & interview_selcomp !=. & cntry == "`x'"  // Cases after (30 days bandwidth)
	local after_bw = `r(N)'
	
	hist interview_selcomp if cntry == "`x'"	,	///
												xaxis(1 2)	///
												freq	///
												xline(`invasion_date', lpattern(dash))	///
												xlabel(, format(%tdDD/NN/YY) axis(1) labsize(vsmall)) 	///
												xlabel(`invasion_date', format(%tdDD/NN/YY) axis(2) notick labsize(vsmall))	///
												xtitle("", axis(2))	///
												xtitle("Date questionnaire finished", axis(1))	///
												ytitle("")	///
												ti("`x'", size(small))
	
	graph save "${aux}interviews_`x'", replace
}

* Generating plot with all the countries
graph combine	"${aux}interviews_Spain"	///
				"${aux}interviews_Poland"	///
				"${aux}interviews_Serbia"	
				
graph export "${output}distribution_interviews_self.eps", replace

**********************************************************
* Variable capturing countries included in the analyses
***********************************************************
gen included = 0

foreach x of local cntries_face{
	replace included = 1 if cntry == "`x'"
}

foreach x of local cntries_self{
	replace included = 1 if cntry == "`x'"
}

tab cntry included, m

replace included = 0 if cntry == "Norway" 

***********************************************************
* Time variable that takes value 0 on day of the invasion
***********************************************************
gen time_zero = . 

foreach x of local cntries_face{
	replace time_zero = interview_face - `invasion_date' if cntry == "`x'"
}

foreach x of local cntries_self{
	replace time_zero = interview_selcomp - `invasion_date' if cntry == "`x'"
}

replace time_zero = . if cntry == "Norway"

****************************************
* Analysis of pre-existing time-trends
****************************************

*** Models to examine trends

* Pooled sample
foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr cttresa votedir gvctzpv grdfinc viepol wpestop keydec{
	reg `var' c.time_zero##c.time_zero##c.time_zero##c.time_zero i.cntrynum if time_zero < 0 & time_zero > -31 
	est store `var'
}

* By country
levelsof cntry if included == 1, local(levels)

foreach x of local levels{
		foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr cttresa votedir  gvctzpv grdfinc viepol wpestop keydec{
		reg `var' c.time_zero##c.time_zero##c.time_zero##c.time_zero if time_zero < 0 & time_zero > -31  & cntry == "`x'"
		est store `var'_`x'
	}
	local `x' = "implvdm_`x' accalaw_`x' stfdem_`x' fairelc_`x'  rghmgpr_`x' votedir_`x' cttresa_`x' gvctzpv_`x' grdfinc_`x' viepol_`x' wpestop_`x' keydec_`x', bylabel(`x')"
}

* Plots to summarize results 
coefplot	(implvdm, label(Importance demo))	///
			(accalaw, label(Strong leader))	///
			(stfdem, label(Satisfaction democracy))	///
			(fairelc, label(Fair elect))	///
			(rghmgpr, label(Minority rights))	///
			(cttresa, label(Rule of law))	///
			(votedir, label(Referenda))	///
			(gvctzpv, label(Poverty))	///
			(grdfinc, label(Income diff))	///
			(viepol, label(Views people))	///
			(wpestop, label(Will people))	///
			(keydec, label(Decisions EU)),	///
			bylabel("Pooled")  ///
			|| `Greece'	///
			|| `Italy'	///
			|| `Macedonia'	///
			|| `Montenegro'	///
			|| `Netherlands'	///
			|| 	, drop(_cons *cntrynum) xline(0, lpattern(solid)) legend(rows(2) size(1.5)) msize(vsmall)	///
				levels(99 95)	///
				coeflabels(time_zero = "Days to invasion"	///
				c.time_zero#c.time_zero = "Days to invasion{superscript: 2}"	///
				c.time_zero#c.time_zero#c.time_zero = "Days to invasion{superscript: 3}"	///
				c.time_zero#c.time_zero#c.time_zero#c.time_zero = "Days to invasion{superscript: 4}")
				
graph export "${output}trends1.eps", replace 

			
coefplot	(implvdm_Poland, label(Importance demo))	///
			(accalaw_Poland, label(Strong leader))	///
			(stfdem_Poland, label(Satisfaction democracy))	///
			(fairelc_Poland, label(Fair elect))	///
			(rghmgpr_Poland, label(Minority rights))	///
			(cttresa_Poland, label(Rule of law))	///
			(votedir_Poland, label(Referenda))	///
			(gvctzpv_Poland, label(Poverty))	///
			(grdfinc_Poland, label(Income diff))	///
			(viepol_Poland, label(Views people))	///
			(wpestop_Poland, label(Will people))	///
			(keydec_Poland, label(Decisions EU)),	///
			bylabel("Poland")  ///
			|| `Portugal'	///
			|| `Serbia'	///
			|| `Spain'	///
			|| `Switzerland'	///
			|| 	, drop(_cons *cntrynum) xline(0, lpattern(solid)) legend(rows(2) size(1.5)) msize(vsmall)	///
				levels(99 95)	///
				coeflabels(time_zero = "Days to invasion"	///
				c.time_zero#c.time_zero = "Days to invasion{superscript: 2}"	///
				c.time_zero#c.time_zero#c.time_zero = "Days to invasion{superscript: 3}"	///
				c.time_zero#c.time_zero#c.time_zero#c.time_zero = "Days to invasion{superscript: 4}")			
			
graph export "${output}trends2.eps", replace 


*** Models that partition sample at mid-point of the control group 
gen placebo_trends = . 
replace placebo_trends = 0 if time_zero < -15 & time_zero > -31 
replace placebo_trends = 1 if time_zero < 0 & time_zero >= -15 


* Pooled sample
foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr cttresa votedir gvctzpv grdfinc viepol wpestop keydec{
	reg `var' placebo_trends i.cntrynum
	est store `var'
}

* By country
levelsof cntry if included == 1, local(levels)

foreach x of local levels{
		foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr cttresa votedir  gvctzpv grdfinc viepol wpestop keydec{
		reg `var' placebo_trends if cntry == "`x'"
		est store `var'_`x'
	}
	local `x' = "implvdm_`x' accalaw_`x' stfdem_`x' fairelc_`x' rghmgpr_`x' votedir_`x' cttresa_`x' gvctzpv_`x' grdfinc_`x' viepol_`x' wpestop_`x' keydec_`x', bylabel(`x')"
}

* Plots to summarize results 
coefplot	(implvdm, label(Importance demo))	///
			(accalaw, label(Strong leader))	///
			(stfdem, label(Satisfaction democracy))	///
			(fairelc, label(Fair elect))	///
			(rghmgpr, label(Minority rights))	///
			(cttresa, label(Rule of law))	///
			(votedir, label(Referenda))	///
			(gvctzpv, label(Poverty))	///
			(grdfinc, label(Income diff))	///
			(viepol, label(Views people))	///
			(wpestop, label(Will people))	///
			(keydec, label(Decisions EU)),	///
			bylabel("Pooled")  ///
			|| `Greece'	///
			|| `Italy'	///
			|| `Macedonia'	///
			|| `Montenegro'	///
			|| `Netherlands'	///
			|| `Poland'	///			
			|| `Portugal'	///
			|| `Serbia'	///
			|| `Spain'	///
			|| `Switzerland'	///
			|| 	, drop(_cons *cntrynum) xline(0, lpattern(solid)) legend(rows(2)) msize(vsmall)	///
				levels(99 95)	///
				coeflabels(placebo_trends = "Placebo cutoff")
				
graph export "${output}placebo1.eps", replace 
				
**************************************************************************
*****++++++++++++++++++++++Treatment definition+++++++++++++++++++++******

*****+++++++++++++++++++++Face-to-face countries+++++++++++++++++++++*****
**************************************************************************

***********************************************************************
* Power analysis by country based on key variable implvdm
***********************************************************************

foreach c of local cntries_face{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_face < `invasion_date' 
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'" 
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'" 
		replace power_treatment = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"  
		replace power_treatment = . if interview_face == `cut' & cntry == "`c'"  
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Dropping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_face{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Switzerland.gph"	///
			"${aux}power_Greece.gph"	///
			"${aux}power_Italy.gph"	///
			"${aux}power_Montenegro.gph"	///
			"${aux}power_Macedonia.gph"	///
			"${aux}power_Netherlands.gph"	///
			"${aux}power_Norway.gph"	///
			"${aux}power_Portugal.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 				
				
graph export "${output}power_face_to_face.eps", replace 

* Returning to master frame
frame change default 
	
**********************************************************
* Treatment variables based on power analyses (by country)
**********************************************************

* Treatment variable for power based on half-standard deviation change 
gen treatment_halfsd = . 
foreach c of local cntries_face{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'"
	replace treatment_halfsd = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"
	replace treatment_halfsd = . if interview_face == `cut' & cntry == "`c'"
	if `bandwidth_half_`c'' == 0{
		replace treatment_halfsd = . if cntry == "`c'"
	}
}

* Treatment variable for power based on one-third-standard deviation change 
gen treatment_thirdsd = . 
foreach c of local cntries_face{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_third_`c''
	local low = (`cut' - `bandwidth_third_`c'') - 1 
	replace treatment_thirdsd = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'"
	replace treatment_thirdsd = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"
	replace treatment_thirdsd = . if interview_face == `cut' & cntry == "`c'"
	if `bandwidth_third_`c'' == 0{
		replace treatment_thirdsd = . if cntry == "`c'"
	}
}


*****************************************************************************
*****++++++++++++++++++++++Treatment definition+++++++++++++++++++++*********

*****+++++++++++++++++++++Self-completion countries+++++++++++++++++++++*****
*****************************************************************************


************* Repeat for the countries that fielded through self-completion 

***********************************************************************
* Power analysis by country based on key variable implvdm
***********************************************************************

foreach c of local cntries_self{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_selcomp < `invasion_date' 
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'" 
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'" 
		replace power_treatment = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"  
		replace power_treatment = . if interview_selcomp == `cut' & cntry == "`c'"  
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Dropping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_self{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Spain.gph"	///
			"${aux}power_Poland.gph"	///
			"${aux}power_Serbia.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 
				
graph export "${output}power_selfcompletion.eps", replace 

* Returning to master frame
frame change default 
	
**********************************************************
* Treatment variables based on power analyses (by country)
**********************************************************

* Treatment variable for power based on half-standard deviation change 
foreach c of local cntries_self{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'"
	replace treatment_halfsd = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"
	replace treatment_halfsd = . if interview_selcomp == `cut' & cntry == "`c'"
	if `bandwidth_half_`c''  == 0{
		replace treatment_halfsd = . if cntry == "`c'"
	}
}

* Treatment variable for power based on one-third-standard deviation change 
foreach c of local cntries_self{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_third_`c''
	local low = (`cut' - `bandwidth_third_`c'') - 1 
	replace treatment_thirdsd = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'"
	replace treatment_thirdsd = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"
	replace treatment_thirdsd = . if interview_selcomp == `cut' & cntry == "`c'"
	if `bandwidth_third_`c'' == 0{
		replace treatment_thirdsd = . if cntry == "`c'"
	}
}


********************
* Minor adjustments 
********************
foreach var of varlist treatment_*{
	replace `var' = . if cntry == "Norway"
}

*****************************************************************************
*****+++++++++++++++++++++++++++++ANALYSIS++++++++++++++++++++++++++*********
*****************************************************************************

**********************************************************
* DESCRIPTION OF TREATMENT AND CONTROL GROUPS BY COUNTRY
**********************************************************

* Size of groups 
foreach var of varlist treatment_*{
	di "NEW TREATMENT VARIBALE `var'"
	bysort cntry: tab `var'
}

* Bandwidth by country and variable
levelsof cntry if treatment_halfsd != ., local(countries)

foreach x of local countries{
	frame change observations_`x'
	sort Bandwidth
	sum Bandwidth if half_sd == 1
	local half_`x' = r(min)
	sum Bandwidth if third_sd == 1
	local third_`x' = r(min)
} 

frame change default

** Locals for plots labels with country name + N + Bandwidth

* Large effects bandwidth 
levelsof cntry if treatment_halfsd != ., local(countries)
foreach x of local countries{
	local band = `half_`x''
	sum treatment_halfsd if cntry == "`x'"
	local observations = r(N)
	local half_lab_`x' = `"`x' = `""`x'" "N = `observations'" "Bandwidth = `band' days""'"'
	local band_plot_`x' = `band'
}

* Small effects bandwidth 
levelsof cntry if treatment_thirdsd != ., local(countries)
foreach x of local countries{
	local band = `third_`x''
	sum treatment_thirdsd if cntry == "`x'"
	local observations = r(N)
	local third_lab_`x' = `"`x' = `""`x'" "N = `observations'" "Bandwidth = `band' days""'"'
}

***************************************
* VARIABLES TO INCLUDE IN THE ANALYSES
***************************************
sum implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

******************************************
* Balance tests: LARGE EFFECTS BANDWIDTH
******************************************

* Reversing treatment variables so that in means comparisons higher numbers means higher (proportion) in treatment group
foreach var of varlist treatment_*{
	gen r_`var' = `var'
	recode r_`var' (1=0) (0=1)
}


** Conducting t-test and storing results in matrix to generate the plot (it is necessary to estimate the confidence intervals "manually") 

levelsof cntry if treatment_halfsd != ., local(countries)

* Generating matrices to store results
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
}

* Estimating t-tests and storing the results for ploting
local i = 0
foreach x of local countries{
	
	local ++ i
	
	if "`x'" == "Spain" | "`x'" == "Poland" | "`x'" == "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted{
		quietly: ttest `var' if cntry == "`x'", by(r_treatment_halfsd) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `var'[1, `i'] = `diff' 
		matrix `var'_CI[1, `i'] = `ll95' 
		matrix `var'_CI[2, `i'] = `ul95' 
		matrix `var'_CI[3, `i'] = `ll90' 
		matrix `var'_CI[4, `i'] = `ul90'
		}
	}
	if "`x'" != "Spain" & "`x'" != "Poland" & "`x'" != "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted  attempts_10 refusals{
			quietly: ttest `var' if cntry == "`x'", by(r_treatment_halfsd) 
			local diff = r(mu_1) - r(mu_2) 
			local degrees = r(df_t)
			local critical_5 = invttail(`degrees', 0.025)
			local confvalue_5 = `critical_5' * r(se)
			local critical_10 = invttail(`degrees', 0.05)
			local confvalue_10 = `critical_10' * r(se)
			local ll95 = `diff' - `confvalue_5'
			local ul95 = `diff' + `confvalue_5'
			local ll90 = `diff' - `confvalue_10'
			local ul90 = `diff' + `confvalue_10'
	
			* Storing main results
			matrix `var'[1, `i'] = `diff' 
			matrix `var'_CI[1, `i'] = `ll95' 
			matrix `var'_CI[2, `i'] = `ul95' 
			matrix `var'_CI[3, `i'] = `ll90' 
			matrix `var'_CI[4, `i'] = `ul90'
		}
	}
}

* Generating plot
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[1] `var'_CI[2]) (`var'_CI[3] `var'_CI[4])))"
}

coefplot	`coefplot_eduyrs_10'	///
			`coefplot_gndr'	///
			`coefplot_age_10'	///
			`coefplot_active'	///
			`coefplot_voted'	///
			`coefplot_attempts_10'	///
			`coefplot_refusals'	///
			,  msize(vsmall)  xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) ///
			xtitle("Mean differences between treatment and control groups")
			
graph export "${output}balance_large_effect.eps", replace	

***** Regions

* Conducting t-tests and generating matrices to store results
foreach x of local countries{
	levelsof region if cntry == "`x'", local(regions)
	
	* Matrices to store results 
	matrix `x' = J(1,`r(r)',.)
	matrix `x'_CI = J(4,`r(r)',.)
	matrix colnames `x' = `regions'
	matrix colnames `x'_CI = `regions'
	
	local c = 0 
	
	foreach r of local regions{
		gen testvar = 1 if region == "`r'"
		replace testvar = 0 if testvar == .
		local ++ c 
		quietly: ttest testvar if cntry == "`x'", by(r_treatment_halfsd) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `x'[1, `c'] = `diff' 
		matrix `x'_CI[1, `c'] = `ll95' 
		matrix `x'_CI[2, `c'] = `ul95' 
		matrix `x'_CI[3, `c'] = `ll90' 
		matrix `x'_CI[4, `c'] = `ul90'
		
		drop testvar	
		
	}
}

* Generating plot syntax 
foreach x of local countries{
	local `x' = "(matrix(`x'), label(`x')  ci((`x'_CI[1] `x'_CI[2]) (`x'_CI[3] `x'_CI[4])))"
}

* Plot 
foreach x of local countries{
	coefplot	``x'' , xline(0, lpattern(solid)) ti("`x'") ylabel(, labsize(vsmall))
	if "`x'" == "Spain"{
		coefplot	``x'' , xline(0, lpattern(solid)) ti("`x'") ylabel(, labsize(tiny))
	}
	graph save "${aux}`x'.gph", replace
}

graph combine 	"${aux}Greece.gph" 	///
				"${aux}Italy.gph"	///
				"${aux}Macedonia.gph"	///
				"${aux}Netherlands.gph"	///
				"${aux}Poland.gph"	///
				"${aux}Portugal.gph"	///
				"${aux}Serbia.gph"	///	
				"${aux}Spain.gph"	///	
				"${aux}Switzerland.gph"	, xcommon
				
graph export "${output}balance_large_effect_regions.eps", replace	
				
******************************************
* Balance tests: SMALL EFFECTS BANDWIDTH
******************************************

** Conducting t-test and storing results in matrix to generate the plot (it is necessary to estimate the confidence intervals "manually") 

levelsof cntry if treatment_thirdsd != ., local(countries)

* Generating matrices to store results
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
}

* Estimating t-tests and storing the results for ploting
local i = 0
foreach x of local countries{
	
	local ++ i
	
	if "`x'" == "Spain" | "`x'" == "Poland" | "`x'" == "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted{
		quietly: ttest `var' if cntry == "`x'", by(treatment_thirdsd) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `var'[1, `i'] = `diff' 
		matrix `var'_CI[1, `i'] = `ll95' 
		matrix `var'_CI[2, `i'] = `ul95' 
		matrix `var'_CI[3, `i'] = `ll90' 
		matrix `var'_CI[4, `i'] = `ul90'
		}
	}
	if "`x'" != "Spain" & "`x'" != "Poland" & "`x'" != "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted  attempts_10 refusals{
		quietly: ttest `var' if cntry == "`x'", by(treatment_thirdsd) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `var'[1, `i'] = `diff' 
		matrix `var'_CI[1, `i'] = `ll95' 
		matrix `var'_CI[2, `i'] = `ul95' 
		matrix `var'_CI[3, `i'] = `ll90' 
		matrix `var'_CI[4, `i'] = `ul90'
		}
	}
}

* Generating plot
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[1] `var'_CI[2]) (`var'_CI[3] `var'_CI[4])))"
}

coefplot	`coefplot_eduyrs_10'	///
			`coefplot_gndr'	///
			`coefplot_age_10'	///
			`coefplot_active'	///
			`coefplot_voted'	///
			`coefplot_attempts_10'	///
			`coefplot_refusals'	///
			,  msize(vsmall)  xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) ///
			xtitle("Mean differences between treatment and control groups")
			
graph export "${output}balance_small_effect.eps", replace	


***** Regions

* Conducting t-tests and generating matrices to store results
foreach x of local countries{
	levelsof region if cntry == "`x'", local(regions)
	
	* Matrices to store results 
	matrix `x' = J(1,`r(r)',.)
	matrix `x'_CI = J(4,`r(r)',.)
	matrix colnames `x' = `regions'
	matrix colnames `x'_CI = `regions'
	
	local c = 0 
	
	foreach r of local regions{
		gen testvar = 1 if region == "`r'"
		replace testvar = 0 if testvar == .
		local ++ c 
		quietly: ttest testvar if cntry == "`x'", by(treatment_thirdsd) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `x'[1, `c'] = `diff' 
		matrix `x'_CI[1, `c'] = `ll95' 
		matrix `x'_CI[2, `c'] = `ul95' 
		matrix `x'_CI[3, `c'] = `ll90' 
		matrix `x'_CI[4, `c'] = `ul90'
		
		drop testvar	
		
	}
}

* Generating plot syntax 
foreach x of local countries{
	local `x' = "(matrix(`x'), label(`x')  ci((`x'_CI[1] `x'_CI[2]) (`x'_CI[3] `x'_CI[4])))"
}

* Plot 
foreach x of local countries{
	coefplot	``x'' , xline(0, lpattern(solid)) ti("`x'") ylabel(, labsize(vsmall))
	graph save "${aux}`x'.gph", replace
}

graph combine 	"${aux}Greece.gph" 	///
				"${aux}Italy.gph"	///
				"${aux}Netherlands.gph"	///
				"${aux}Poland.gph"	///
				"${aux}Serbia.gph"	///	
				"${aux}Spain.gph"	///	
				
graph export "${output}balance_small_effect_regions.eps", replace	

************************************
* Analysis: LARGE EFFECTS BANDWIDTH
************************************
* Matrices to store results 
levelsof cntry if treatment_halfsd != ., local(countries)
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
	
	matrix `var'_corr = J(1,10,.)
	matrix colnames `var'_corr = `countries'
}

* Analysis
local i = 0
foreach x of local countries{
	
	local ++ i 
	
	* Balancing weight definition
	if "`x'" != "Montenegro"{
	ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							if cntry == "`x'",  generate(balance) targets(3) tolerance(.050)
	
	svyset [pweight=balance]
	}
	* Analysis 
	foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'", level(95)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'", level(95)
		}		
		local b = r(table)["b", 1]
		local ll = r(table)["ll", 1]
		local ul = r(table)["ul", 1]
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'", level(99)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'", level(99)
		}	
		local ll99 = r(table)["ll", 1]
		local ul99 = r(table)["ul", 1]
		
	* Storing main results
	matrix `var'[1, `i'] = `b' 
	matrix `var'_CI[1, `i'] = `ll' 
	matrix `var'_CI[2, `i'] = `ul' 
	matrix `var'_CI[3, `i'] = `ll99' 
	matrix `var'_CI[4, `i'] = `ul99'
	}
	
	* Multiple hypothesis testing corrections 
	if "`x'" == "Greece"{
		local num = 9
	} 
	if "`x'" == "Italy"{
		local num = 12
	} 
	if "`x'" == "Macedonia"{
		local num = 14
	} 
	if "`x'" == "Montenegro"{
		local num = 15
	} 
	if "`x'" == "Netherlands"{
		local num = 16
	} 
	if "`x'" == "Poland"{
		local num = 18
	} 
	if "`x'" == "Portugal"{
		local num = 19
	} 
	if "`x'" == "Serbia"{
		local num = 20
	} 
	if "`x'" == "Spain"{
		local num = 23
	} 
	if "`x'" == "Switzerland"{
		local num = 25
	} 	
	
	if "`x'" != "Montenegro"{
		rwolf2 	(svy: reg implvdm treatment_halfsd  if cntrynum == `num')	///
				(svy: reg accalaw treatment_halfsd  if cntrynum == `num') ///
				(svy: reg stfdem treatment_halfsd  if cntrynum == `num')	///
				(svy: reg fairelc treatment_halfsd  if cntrynum == `num') ///
				(svy: reg rghmgpr treatment_halfsd  if cntrynum == `num')	///
				(svy: reg votedir treatment_halfsd  if cntrynum == `num')	///
				(svy: reg cttresa treatment_halfsd  if cntrynum == `num')	///
				(svy: reg gvctzpv treatment_halfsd  if cntrynum == `num')	///
				(svy: reg grdfinc treatment_halfsd  if cntrynum == `num')	///
				(svy: reg viepol treatment_halfsd  if cntrynum == `num')	///
				(svy: reg wpestop treatment_halfsd  if cntrynum == `num')	///
				(svy: reg keydec treatment_halfsd  if cntrynum == `num')	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
		if "`x'" == "Montenegro"{
		rwolf2 	(reg implvdm treatment_halfsd if cntrynum == `num')	///
				(reg accalaw treatment_halfsd if cntrynum == `num') ///
				(reg stfdem treatment_halfsd if cntrynum == `num')	///
				(reg fairelc treatment_halfsd if cntrynum == `num') ///
				(reg rghmgpr treatment_halfsd if cntrynum == `num')	///
				(reg votedir treatment_halfsd if cntrynum == `num')	///
				(reg cttresa treatment_halfsd if cntrynum == `num')	///
				(reg gvctzpv treatment_halfsd if cntrynum == `num')	///
				(reg grdfinc treatment_halfsd if cntrynum == `num')	///
				(reg viepol  treatment_halfsd if cntrynum == `num')	///
				(reg wpestop treatment_halfsd if cntrynum == `num')	///
				(reg keydec  treatment_halfsd if cntrynum == `num')	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
	* Storing resuls multiple hypotheses corrections
	matrix multiple =  e(RW) 
	local r = 0 
	foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		local ++ r 
		matrix `var'_corr[1, `i'] = multiple[`r',3] 
	}
	
	
	* Drop balancing weights variable
	if "`x'" != "Montenegro"{
			drop balance
		}
}

* Plots main results 
foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[3] `var'_CI[4]) (`var'_CI[1] `var'_CI[2])))"
}

coefplot	`coefplot_implvdm'	///
			`coefplot_accalaw'	///
			`coefplot_stfdem'	///
			`coefplot_fairelc'	///
			`coefplot_rghmgpr'	///
			`coefplot_cttresa'	///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			xlabel(-4(1)4)	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')
			
graph export "${output}large_effects_var1.eps", replace
			
coefplot	`coefplot_votedir' ///
			`coefplot_gvctzpv'	///
			`coefplot_grdfinc'	///
			`coefplot_viepol'	///
			`coefplot_wpestop'	///
			`coefplot_keydec'	///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) ///
			xlabel(-4(1)4)	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')
			
graph export "${output}large_effects_var2.eps", replace


* Plots multiple hypothesis testing 
foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	matrix `var' = `var'_corr
	local coefplot_`var' = "(matrix(`var'), label(`label'))"
	matrix colnames `var' = `countries'
}

coefplot	`coefplot_implvdm'	///
			`coefplot_accalaw'	///
			`coefplot_stfdem'	///
			`coefplot_fairelc'	///
			`coefplot_rghmgpr'	///
			`coefplot_cttresa'	///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
			xtitle(" Romano-Wolf stepdown p-values for multiple hypothesis testing")	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')
			
graph export "${output}multiple_correction_var1.eps", replace


coefplot	`coefplot_votedir' ///
			`coefplot_gvctzpv'	///
			`coefplot_grdfinc'	///
			`coefplot_viepol'	///
			`coefplot_wpestop'	///
			`coefplot_keydec'	///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
			xtitle(" Romano-Wolf stepdown p-values for multiple hypothesis testing")	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')
			
graph export "${output}multiple_correction_var2.eps", replace
		
************************************
* Analysis: SMALL EFFECTS BANDWIDTH
************************************
matrix drop _all

* Matrices to store results 
levelsof cntry if treatment_thirdsd != ., local(countries)
foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
	
	matrix `var'_corr = J(1,10,.)
	matrix colnames `var'_corr = `countries'
}

* Analysis
local i = 0
foreach x of local countries{
	
	local ++ i 
	
	* Balancing weight
	if "`x'" != "Montenegro"{
	ebalance treatment_thirdsd	eduyrs gndr agea active	voted						///
							if cntry == "`x'",  generate(balance) targets(3) tolerance(.050)
	
	svyset [pweight=balance]
	}
	* Analysis 
	foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_thirdsd  if cntry == "`x'", level(95)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_thirdsd if cntry == "`x'", level(95)
		}		
		local b = r(table)["b", 1]
		local ll = r(table)["ll", 1]
		local ul = r(table)["ul", 1]
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_thirdsd  if cntry == "`x'", level(99)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_thirdsd if cntry == "`x'", level(99)
		}	
		local ll99 = r(table)["ll", 1]
		local ul99 = r(table)["ul", 1]
		
	* Storing results
	matrix `var'[1, `i'] = `b' 
	matrix `var'_CI[1, `i'] = `ll' 
	matrix `var'_CI[2, `i'] = `ul' 
	matrix `var'_CI[3, `i'] = `ll99' 
	matrix `var'_CI[4, `i'] = `ul99'
	}
	
	* Multiple hypothesis testing corrections 
	if "`x'" == "Greece"{
		local num = 9
	} 
	if "`x'" == "Italy"{
		local num = 12
	} 
	if "`x'" == "Macedonia"{
		local num = 14
	} 
	if "`x'" == "Montenegro"{
		local num = 15
	} 
	if "`x'" == "Netherlands"{
		local num = 16
	} 
	if "`x'" == "Poland"{
		local num = 18
	} 
	if "`x'" == "Portugal"{
		local num = 19
	} 
	if "`x'" == "Serbia"{
		local num = 20
	} 
	if "`x'" == "Spain"{
		local num = 23
	} 
	if "`x'" == "Switzerland"{
		local num = 25
	} 	
	
	if "`x'" != "Montenegro"{
		rwolf2 	(svy: reg implvdm treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg accalaw treatment_thirdsd  if cntrynum == `num') ///
				(svy: reg stfdem treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg fairelc treatment_thirdsd  if cntrynum == `num') ///
				(svy: reg rghmgpr treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg votedir treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg cttresa treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg gvctzpv treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg grdfinc treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg viepol treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg wpestop treatment_thirdsd  if cntrynum == `num')	///
				(svy: reg keydec treatment_thirdsd  if cntrynum == `num')	///
				, indepvars(treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
		if "`x'" == "Montenegro"{
		rwolf2 	(reg implvdm  treatment_thirdsd if cntrynum == `num')	///
				(reg accalaw  treatment_thirdsd if cntrynum == `num') ///
				(reg stfdem  treatment_thirdsd if cntrynum == `num')	///
				(reg fairelc  treatment_thirdsd if cntrynum == `num') ///
				(reg rghmgpr  treatment_thirdsd if cntrynum == `num')	///
				(reg votedir  treatment_thirdsd if cntrynum == `num')	///
				(reg cttresa  treatment_thirdsd if cntrynum == `num')	///
				(reg gvctzpv  treatment_thirdsd if cntrynum == `num')	///
				(reg grdfinc  treatment_thirdsd if cntrynum == `num')	///
				(reg viepol  treatment_thirdsd  if cntrynum == `num')	///
				(reg wpestop  treatment_thirdsd if cntrynum == `num')	///
				(reg keydec  treatment_thirdsd  if cntrynum == `num')	///
				, indepvars(treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
	* Storing resuls multiple hypotheses corrections
	matrix multiple =  e(RW) 
	local r = 0 
	foreach var of varlist implvdm accalaw stfdem fairelc rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		local ++ r 
		matrix `var'_corr[1, `i'] = multiple[`r',3] 
	}
	
	
	* Drop balancing weights variable
	if "`x'" != "Montenegro"{
			drop balance
	}
}

* Combined plot 
* Regular plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[3] `var'_CI[4]) (`var'_CI[1] `var'_CI[2])))"
}

* Correction plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	matrix `var'_corr = `var'_corr
	local coefplot_`var'_corr = "(matrix(`var'_corr), label(`label'))"
	matrix colnames `var' = `countries'
}

* Plot for first group of variables
coefplot	`coefplot_implvdm' `coefplot_accalaw' `coefplot_stfdem' `coefplot_fairelc' `coefplot_rghmgpr' `coefplot_cttresa' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_implvdm_corr' `coefplot_accalaw_corr' `coefplot_stfdem_corr' `coefplot_fairelc_corr' 	`coefplot_rghmgpr_corr' `coefplot_cttresa_corr'	///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law"))	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))
				
graph export "${output}small_effects_var1.eps", replace 
			
* Plot for second group of variables
coefplot	`coefplot_votedir'  `coefplot_gvctzpv' `coefplot_grdfinc' `coefplot_viepol'	`coefplot_wpestop' `coefplot_keydec' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_votedir_corr'  `coefplot_gvctzpv_corr' `coefplot_grdfinc_corr' `coefplot_viepol_corr'	`coefplot_wpestop_corr' `coefplot_keydec_corr' ///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))

graph export "${output}small_effects_var2.eps"	, replace 		
				

*****************************************
* Analysis: POOLED SAMPLE LARGE EFFECTS
*****************************************

* Matrices to store results 

matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balacing weight 
ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis main effects
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' treatment_halfsd i.cntrynum, level(95)
	local b = r(table)["b", 1]
	local ll = r(table)["ll", 1]
	local ul = r(table)["ul", 1]
	svy: reg `var' treatment_halfsd i.cntrynum, level(99)
	local ll99 = r(table)["ll", 1]
	local ul99 = r(table)["ul", 1]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}


* Analsysis multiple hypothesis test corrections
rwolf2 	(svy: reg implvdm treatment_halfsd i.cntrynum)	///
		(svy: reg accalaw treatment_halfsd i.cntrynum) ///
		(svy: reg stfdem treatment_halfsd i.cntrynum)	///
		(svy: reg fairelc treatment_halfsd i.cntrynum) ///
		(svy: reg rghmgpr treatment_halfsd i.cntrynum)	///
		(svy: reg votedir treatment_halfsd i.cntrynum)	///
		(svy: reg cttresa treatment_halfsd i.cntrynum)	///
		(svy: reg gvctzpv treatment_halfsd i.cntrynum)	///
		(svy: reg grdfinc treatment_halfsd i.cntrynum)	///
		(svy: reg viepol treatment_halfsd i.cntrynum)	///
		(svy: reg wpestop treatment_halfsd i.cntrynum)	///
		(svy: reg keydec treatment_halfsd i.cntrynum)	///
		, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
		reps(`replications') verbose seed(`seed_corrections')

matrix multiple =  e(RW) 		
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	matrix results_corr[1,`c'] = multiple[`c',3] 
}


drop balance


* Plot main results
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	

graph export "${output}large_effects_pooled.eps", replace


* Plot p-values for multiple hypothesis testing correction
coefplot	(matrix(results_corr)) ///
			,  msize(small) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
			xtitle(" Romano-Wolf stepdown p-values for multiple hypothesis testing")	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	

graph export "${output}correction_pooled.eps", replace

*******************************************************
* Analysis of potential interaction with news exposure 
*******************************************************

*** Data management (dichotomous measure of political interest)
gen news_exposure = . 

gen news_exposure2 = . 

sum nwspol if included == 1, det

replace news_exposure = 0 if nwspol == 0 & nwspol != .

replace news_exposure = 1 if nwspol > 0 & nwspol != .

replace news_exposure2 = 0 if nwspol <= `r(p50)' & nwspol != .

replace news_exposure2 = 1 if nwspol > `r(p50)' & nwspol != .

*** Pooled sample

* Matrices to store results 

matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balancing weight 
ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis main effects (news exposure (Yes vs No))
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' i.treatment_halfsd##i.news_exposure i.cntrynum, level(95)
	local b = r(table)["b", 8]
	local ll = r(table)["ll", 8]
	local ul = r(table)["ul", 8]
	svy: reg `var' treatment_halfsd##i.news_exposure i.cntrynum, level(99)
	local ll99 = r(table)["ll", 8]
	local ul99 = r(table)["ul", 8]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}

* Plot main results
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)
			
graph export "${output}interaction_news_exposure_yes_no.eps", replace

* Analysis main effects (news exposure (Lower or higher than median))
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' i.treatment_halfsd##i.news_exposure2 i.cntrynum, level(95)
	local b = r(table)["b", 8]
	local ll = r(table)["ll", 8]
	local ul = r(table)["ul", 8]
	svy: reg `var' treatment_halfsd##i.news_exposure2 i.cntrynum, level(99)
	local ll99 = r(table)["ll", 8]
	local ul99 = r(table)["ul", 8]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}

* Plot main results
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			,  msize(vsmall)  ciopts(lwidth(thin medium)) xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)
			
graph export "${output}interaction_news_exposure_median.eps", replace



drop balance


***************************************
* Analysis of non-response (atrittion)
***************************************
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	gen miss_`var' = `var' 
	recode miss_`var' (0/10=0) (.a=1) (.b=1) (.c=1)
}

* Pooled sample
foreach var of varlist miss_*{
	reg `var' treatment_halfsd i.cntrynum
	est store `var'
}

* By country
levelsof cntry if included == 1, local(levels)

foreach x of local levels{
		foreach var of varlist miss_*{
		reg `var' treatment_halfsd if cntry == "`x'"
		est store `var'_`x'
	}
	local `x' = "miss_implvdm_`x' miss_accalaw_`x' miss_stfdem_`x' miss_fairelc_`x'  miss_rghmgpr_`x' miss_votedir_`x' miss_cttresa_`x' miss_gvctzpv_`x' miss_grdfinc_`x' miss_viepol_`x' miss_wpestop_`x' miss_keydec_`x', bylabel(`x')"
}

* Plots to summarize results 
coefplot	(miss_implvdm, label(Importance demo))	///
			(miss_accalaw, label(Strong leader))	///
			(miss_stfdem, label(Satisfaction democracy))	///
			(miss_fairelc, label(Fair elect))	///
			(miss_rghmgpr, label(Minority rights))	///
			(miss_cttresa, label(Rule of law))	///
			(miss_votedir, label(Referenda))	///
			(miss_gvctzpv, label(Poverty))	///
			(miss_grdfinc, label(Income diff))	///
			(miss_viepol, label(Views people))	///
			(miss_wpestop, label(Will people))	///
			(miss_keydec, label(Decisions EU)),	///
			bylabel("Pooled")  ///
			|| `Greece'	///
			|| `Italy'	///
			|| `Macedonia'	///
			|| `Montenegro'	///
			|| `Netherlands'	///
			|| `Poland'	///			
			|| `Portugal'	///
			|| `Serbia'	///
			|| `Spain'	///
			|| `Switzerland'	///
			|| 	, drop(_cons *cntrynum) xline(0, lpattern(solid)) legend(rows(2)) msize(vsmall)	///
				levels(99 95)	///
				coeflabels(treatment_halfsd = "Treatment effect")
				
graph export "${output}missing.eps", replace 

*****************************************
* Analysis: POOLED SAMPLE SMALL EFFECTS
*****************************************
* Matrices to store results 
matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balacing weight 
ebalance treatment_thirdsd	eduyrs gndr agea active	voted						///
							,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis 
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' treatment_thirdsd i.cntrynum, level(95)
	local b = r(table)["b", 1]
	local ll = r(table)["ll", 1]
	local ul = r(table)["ul", 1]
	svy: reg `var' treatment_thirdsd i.cntrynum, level(99)
	local ll99 = r(table)["ll", 1]
	local ul99 = r(table)["ul", 1]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}

* Analsysis multiple hypothesis test corrections
rwolf2 	(svy: reg implvdm treatment_thirdsd i.cntrynum)	///
		(svy: reg accalaw treatment_thirdsd i.cntrynum) ///
		(svy: reg stfdem treatment_thirdsd i.cntrynum)	///
		(svy: reg fairelc treatment_thirdsd i.cntrynum) ///
		(svy: reg rghmgpr treatment_thirdsd i.cntrynum)	///
		(svy: reg votedir treatment_thirdsd i.cntrynum)	///
		(svy: reg cttresa treatment_thirdsd i.cntrynum)	///
		(svy: reg gvctzpv treatment_thirdsd i.cntrynum)	///
		(svy: reg grdfinc treatment_thirdsd i.cntrynum)	///
		(svy: reg viepol treatment_thirdsd i.cntrynum)	///
		(svy: reg wpestop treatment_thirdsd i.cntrynum)	///
		(svy: reg keydec treatment_thirdsd i.cntrynum)	///
		, indepvars(treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd, treatment_thirdsd)	///
		reps(`replications') verbose seed(`seed_corrections')

matrix multiple =  e(RW) 		
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	matrix results_corr[1,`c'] = multiple[`c',3] 
}



drop balance

* Combined plot 
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			, bylabel((A) Estimation)	///
			|| (matrix(results_corr)) ///
			, bylabel((B) Multiple testing corrections)	///
			||,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	///
			byopts(xrescale)

addplot 1: ,	legend(off) 	///
				xlabel(-1(1)1) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(off) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))


graph export "${output}small_effects_pooled.eps", replace


*************************************************
* Analysis: OVERPOWERED SAMPLE IN ALL COUNTRIES
************************************************
matrix drop _all

* Treatment definition based on the same bandwidth (9 days) in each countrys

gen treatment_power = . 

* Face-to-face
foreach c of local cntries_face{
	local cut = `invasion_date'
	local up = `cut' +  30
	local low = (`cut' - 30) - 1 
	replace treatment_power = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'"
	replace treatment_power = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"
	replace treatment_power = . if interview_face == `cut' & cntry == "`c'"
}

* Self-completion
foreach c of local cntries_self{
	local cut = `invasion_date'
	local up = `cut' + 30
	local low = (`cut' - 30) - 1 
	replace treatment_power = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'"
	replace treatment_power = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"
	replace treatment_power = . if interview_selcomp == `cut' & cntry == "`c'"
}


* We drop Norway for comparability (and lack of power)
foreach var of varlist treatment_*{
	replace `var' = . if cntry == "Norway"
}

** Locals for plots labels with country name + N + Bandwidth

levelsof cntry if treatment_power != ., local(countries)
foreach x of local countries{
	local band =  30
	sum treatment_power if cntry == "`x'"
	local observations = r(N)
	local half_lab_`x' = `"`x' = `""`x'" "N = `observations'" "Bandwidth = `band' days""'"'
	local band_plot_`x' = `band'
}

* Matrices to store results 
levelsof cntry if treatment_power != ., local(countries)
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
	
	matrix `var'_corr = J(1,10,.)
	matrix colnames `var'_corr = `countries'
}

* Analysis
local i = 0
foreach x of local countries{
	
	local ++ i 
	
	* Balancing weight
	if "`x'" != "Montenegro"{
	ebalance treatment_power	eduyrs gndr agea active	voted						///
							if cntry == "`x'",  generate(balance) targets(3) tolerance(.050)
	
	svyset [pweight=balance]
	}
	* Analysis 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_power  if cntry == "`x'", level(95)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_power if cntry == "`x'", level(95)
		}		
		local b = r(table)["b", 1]
		local ll = r(table)["ll", 1]
		local ul = r(table)["ul", 1]
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_power  if cntry == "`x'", level(99)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_power if cntry == "`x'", level(99)
		}	
		local ll99 = r(table)["ll", 1]
		local ul99 = r(table)["ul", 1]
		
	* Storing results
	matrix `var'[1, `i'] = `b' 
	matrix `var'_CI[1, `i'] = `ll' 
	matrix `var'_CI[2, `i'] = `ul' 
	matrix `var'_CI[3, `i'] = `ll99' 
	matrix `var'_CI[4, `i'] = `ul99'
	}
	
	* Multiple hypothesis testing corrections 
	if "`x'" == "Greece"{
		local num = 9
	} 
	if "`x'" == "Italy"{
		local num = 12
	} 
	if "`x'" == "Macedonia"{
		local num = 14
	} 
	if "`x'" == "Montenegro"{
		local num = 15
	} 
	if "`x'" == "Netherlands"{
		local num = 16
	} 
	if "`x'" == "Poland"{
		local num = 18
	} 
	if "`x'" == "Portugal"{
		local num = 19
	} 
	if "`x'" == "Serbia"{
		local num = 20
	} 
	if "`x'" == "Spain"{
		local num = 23
	} 
	if "`x'" == "Switzerland"{
		local num = 25
	} 	
	
	if "`x'" != "Montenegro"{
		rwolf2 	(svy: reg implvdm treatment_power  if cntrynum == `num')	///
				(svy: reg accalaw treatment_power  if cntrynum == `num') ///
				(svy: reg stfdem treatment_power  if cntrynum == `num')	///
				(svy: reg fairelc treatment_power  if cntrynum == `num') ///
				(svy: reg rghmgpr treatment_power  if cntrynum == `num')	///
				(svy: reg votedir treatment_power  if cntrynum == `num')	///
				(svy: reg cttresa treatment_power  if cntrynum == `num')	///
				(svy: reg gvctzpv treatment_power  if cntrynum == `num')	///
				(svy: reg grdfinc treatment_power  if cntrynum == `num')	///
				(svy: reg viepol treatment_power  if cntrynum == `num')	///
				(svy: reg wpestop treatment_power  if cntrynum == `num')	///
				(svy: reg keydec treatment_power  if cntrynum == `num')	///
				, indepvars(treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
		if "`x'" == "Montenegro"{
		rwolf2 	(reg implvdm  treatment_power if cntrynum == `num')	///
				(reg accalaw  treatment_power if cntrynum == `num') ///
				(reg stfdem  treatment_power if cntrynum == `num')	///
				(reg fairelc  treatment_power if cntrynum == `num') ///
				(reg rghmgpr  treatment_power if cntrynum == `num')	///
				(reg votedir  treatment_power if cntrynum == `num')	///
				(reg cttresa  treatment_power if cntrynum == `num')	///
				(reg gvctzpv  treatment_power if cntrynum == `num')	///
				(reg grdfinc  treatment_power if cntrynum == `num')	///
				(reg viepol  treatment_power  if cntrynum == `num')	///
				(reg wpestop  treatment_power if cntrynum == `num')	///
				(reg keydec  treatment_power  if cntrynum == `num')	///
				, indepvars(treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
	* Storing resuls multiple hypotheses corrections
	matrix multiple =  e(RW) 
	local r = 0 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		local ++ r 
		matrix `var'_corr[1, `i'] = multiple[`r',3] 
	}
	
	
	* Drop balancing weights variable
	if "`x'" != "Montenegro"{
			drop balance
	}
}

* Combined plot 
* Regular plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[3] `var'_CI[4]) (`var'_CI[1] `var'_CI[2])))"
}

* Corrections plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	matrix `var'_corr = `var'_corr
	local coefplot_`var'_corr = "(matrix(`var'_corr), label(`label'))"
	matrix colnames `var' = `countries'
}

* Plot for first group of variables
coefplot	`coefplot_implvdm' `coefplot_accalaw' `coefplot_stfdem' `coefplot_fairelc' 	`coefplot_rghmgpr' `coefplot_cttresa' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_implvdm_corr' `coefplot_accalaw_corr' `coefplot_stfdem_corr' `coefplot_fairelc_corr' 	`coefplot_rghmgpr_corr' `coefplot_cttresa_corr'	///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law"))	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))
				
graph export "${output}treatment_power_var1.eps", replace
			
* Plot for second group of variables
coefplot	`coefplot_votedir'  `coefplot_gvctzpv' `coefplot_grdfinc' `coefplot_viepol'	`coefplot_wpestop' `coefplot_keydec' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_votedir_corr'  `coefplot_gvctzpv_corr' `coefplot_grdfinc_corr' `coefplot_viepol_corr'	`coefplot_wpestop_corr' `coefplot_keydec_corr' ///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`half_lab_Greece' `half_lab_Italy' `half_lab_Macedonia' `half_lab_Montenegro' `half_lab_Netherlands' `half_lab_Poland' `half_lab_Portugal' `half_lab_Serbia' `half_lab_Spain' `half_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))

graph export "${output}treatment_power_var2.eps", replace
				
****** Pooled sample 
* Matrices to store results 

matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balacing weight 
ebalance treatment_power	eduyrs gndr agea active	voted						///
							,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis main effects
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' treatment_power i.cntrynum, level(95)
	local b = r(table)["b", 1]
	local ll = r(table)["ll", 1]
	local ul = r(table)["ul", 1]
	svy: reg `var' treatment_power i.cntrynum, level(99)
	local ll99 = r(table)["ll", 1]
	local ul99 = r(table)["ul", 1]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}


*Analsysis multiple hypothesis test corrections
rwolf2 	(svy: reg implvdm treatment_power i.cntrynum)	///
		(svy: reg accalaw treatment_power i.cntrynum) ///
		(svy: reg stfdem treatment_power i.cntrynum)	///
		(svy: reg fairelc treatment_power i.cntrynum) ///
		(svy: reg rghmgpr treatment_power i.cntrynum)	///
		(svy: reg votedir treatment_power i.cntrynum)	///
		(svy: reg cttresa treatment_power i.cntrynum)	///
		(svy: reg gvctzpv treatment_power i.cntrynum)	///
		(svy: reg grdfinc treatment_power i.cntrynum)	///
		(svy: reg viepol treatment_power i.cntrynum)	///
		(svy: reg wpestop treatment_power i.cntrynum)	///
		(svy: reg keydec treatment_power i.cntrynum)	///
		, indepvars(treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power, treatment_power)	///
		reps(`replications') verbose seed(`seed_corrections')

matrix multiple =  e(RW) 		
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	matrix results_corr[1,`c'] = multiple[`c',3] 
}


drop balance


* Combined plot 
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			, bylabel((A) Estimation)	///
			|| (matrix(results_corr)) ///
			, bylabel((B) Multiple testing corrections)	///
			||,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	///
			byopts(xrescale)

addplot 1: ,	legend(off) 	///
				xlabel(-1(1)1) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(off) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))


graph export "${output}treatment_power_pooled.eps", replace

*** Balance tests for overpowered sample 

* Reversing treatment variables so that in means comparisons higher numbers means higher (proportion) in treatment group
foreach var of varlist treatment_power{
	gen r_`var' = `var'
	recode r_`var' (1=0) (0=1)
}


** Conducting t-test and storing results in matrix to generate the plot (it is necessary to estimate the confidence intervals "manually") 

levelsof cntry if treatment_power != ., local(countries)

* Generating matrices to store results
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
}


* Estimating t-tests and storing the results for ploting
local i = 0
foreach x of local countries{
	
	local ++ i
	
	if "`x'" == "Spain" | "`x'" == "Poland" | "`x'" == "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted{
		quietly: ttest `var' if cntry == "`x'", by(treatment_power) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `var'[1, `i'] = `diff' 
		matrix `var'_CI[1, `i'] = `ll95' 
		matrix `var'_CI[2, `i'] = `ul95' 
		matrix `var'_CI[3, `i'] = `ll90' 
		matrix `var'_CI[4, `i'] = `ul90'
		}
	}
	if "`x'" != "Spain" & "`x'" != "Poland" & "`x'" != "Serbia"{
		foreach var of varlist eduyrs_10 gndr age_10 active voted  attempts_10 refusals{
		quietly: ttest `var' if cntry == "`x'", by(treatment_power) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `var'[1, `i'] = `diff' 
		matrix `var'_CI[1, `i'] = `ll95' 
		matrix `var'_CI[2, `i'] = `ul95' 
		matrix `var'_CI[3, `i'] = `ll90' 
		matrix `var'_CI[4, `i'] = `ul90'
		}
	}
}

* Generating plot
foreach var of varlist eduyrs_10 gndr age_10 active voted attempts_10 refusals{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[1] `var'_CI[2]) (`var'_CI[3] `var'_CI[4])))"
}

coefplot	`coefplot_eduyrs_10'	///
			`coefplot_gndr'	///
			`coefplot_age_10'	///
			`coefplot_active'	///
			`coefplot_voted'	///
			`coefplot_attempts_10'	///
			`coefplot_refusals'	///
			,  msize(vsmall)  xline(0, lpattern(solid)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) ///
			xtitle("Mean differences between treatment and control groups")
			
graph export "${output}balance_overpowered_sample.eps", replace	


***** Regions

* Conducting t-tests and generating matrices to store results
foreach x of local countries{
	levelsof region if cntry == "`x'", local(regions)
	
	* Matrices to store results 
	matrix `x' = J(1,`r(r)',.)
	matrix `x'_CI = J(4,`r(r)',.)
	matrix colnames `x' = `regions'
	matrix colnames `x'_CI = `regions'
	
	local c = 0 
	
	foreach r of local regions{
		gen testvar = 1 if region == "`r'"
		replace testvar = 0 if testvar == .
		local ++ c 
		quietly: ttest testvar if cntry == "`x'", by(treatment_power) 
		local diff = r(mu_1) - r(mu_2) 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
	
		* Storing main results
		matrix `x'[1, `c'] = `diff' 
		matrix `x'_CI[1, `c'] = `ll95' 
		matrix `x'_CI[2, `c'] = `ul95' 
		matrix `x'_CI[3, `c'] = `ll90' 
		matrix `x'_CI[4, `c'] = `ul90'
		
		drop testvar	
		
	}
}

* Generating plot syntax 
foreach x of local countries{
	local `x' = "(matrix(`x'), label(`x')  ci((`x'_CI[1] `x'_CI[2]) (`x'_CI[3] `x'_CI[4])))"
}

* Plot 
foreach x of local countries{
	coefplot	``x'' , xline(0, lpattern(solid)) ti("`x'") ylabel(, labsize(vsmall))
	if "`x'" == "Spain"{
		coefplot	``x'' , xline(0, lpattern(solid)) ti("`x'") ylabel(, labsize(tiny))
	}
	graph save "${aux}`x'.gph", replace
}

graph combine 	"${aux}Greece.gph" 	///
				"${aux}Italy.gph"	///
				"${aux}Macedonia.gph"	///
				"${aux}Netherlands.gph"	///
				"${aux}Poland.gph"	///
				"${aux}Portugal.gph"	///
				"${aux}Serbia.gph"	///	
				"${aux}Spain.gph"	///	
				"${aux}Switzerland.gph"	, xcommon
				
graph export "${output}balance_treatment_power_regions.eps", replace	


******************************************
*** Salience of the event. Google Trends
******************************************

*** Worldwide searches 
frame create google_world 
frame change google_world

import delimited "${gtrends}world.csv", delimiter(comma) clear 

drop in 1 
drop in 1

rename v2 searches

replace searches = "0" if searches == "<1"
destring searches, replace

generate date = date(v1,"YMD")
label var date "Date"
format date %tddd/nn/YY

local jan = date("January 1, 2022","MDY")
local feb = date("February 1, 2022","MDY")
local mar = date("March 31, 2022","MDY")
loca cut = date("February 24, 2022","MDY")

graph twoway 	(line searches date,	///
				xline(`cut', lpattern(solid))	///
				xtitle("Date") xlabel(`jan'  `feb' `cut' `mar')	///
				ytitle("")	///
				ti("Worldwide"))
				
graph save "${aux}gtrends_world.gph", replace

*** Country specific searches
frame change default
levelsof cntry if treatment_halfsd != ., local(countries)
foreach x of local countries{
	di "`x'"
	frame create `x'
	frame change `x'
	import delimited "${gtrends}`x'.csv", delimiter(comma) clear 
	drop in 1 
	drop in 1

	rename v2 searches
	rename v3 pm

	label var searches "Relative daily Google searches for in period (max=100)"

	replace searches = "0" if searches == "<1"
	replace pm = "0" if pm == "<1"

	destring searches, replace
	destring pm, replace

	generate date = date(v1,"YMD")
	label var date "Date"
	format date %tddd/nn/YY

	local jan = date("January 1, 2022","MDY")
	local feb = date("February 1, 2022","MDY")
	local mar = date("March 31, 2022","MDY")
	loca cut = date("February 24, 2022","MDY")
	local ll = date("February 24, 2022","MDY") - `band_plot_`x''
	local ul = date("February 24, 2022","MDY") + `band_plot_`x''
	
	graph twoway 	(line searches date,	///
					xline(`cut', lpattern(solid))	///
					xline(`ll', lpattern(dash) lcolor(green))	///
					xline(`ul', lpattern(dash) lcolor(green))	///
					xtitle("Date")	///
					xlabel(`jan'  `feb' `cut' `mar')) ///
					(line pm date, ///
					legend(ti("{bf:Search terms}") order(1 "Ukraine" 2 "Prime Minister of country") rows(1) pos(6)) ///
					ti("`x'"))

	graph save "${aux}gtrends_`x'.gph", replace 
	frame change default
	frame drop `x'
}

grc1leg2 		"${aux}gtrends_world.gph"	///
				"${aux}gtrends_Greece.gph"	///
				"${aux}gtrends_Italy.gph"	///
				"${aux}gtrends_Macedonia.gph"	///
				"${aux}gtrends_Montenegro.gph"	///
				"${aux}gtrends_Netherlands.gph"	///
				"${aux}gtrends_Poland.gph"	///
				"${aux}gtrends_Portugal.gph"	///
				"${aux}gtrends_Serbia.gph"	///
				"${aux}gtrends_Spain.gph"	///
				"${aux}gtrends_Switzerland.gph"	///
				, ycommon			///
				l1("Relative daily Google searches per country in period (max = 100)", size(small)) ///
				legendfrom("${aux}gtrends_Greece.gph")  
				
	
*graph save "${output}study2_gtrends_countries.gph", replace
graph export "${output}gtrends.eps", replace
								
*****************************************
* Alternative treatment date Februrary 22
*****************************************
* Defining new date
matrix drop _all

levelsof cntry if treatment_power != ., local(countries)
foreach x of local countries{
	frame drop observations_`x'
}

frame drop observations_Norway

drop treatment_*

* Locals for date and initial bandwidth arround invasion date (1 month)
local invasion_date = date("22feb2022","DMY")
local invasion_ll = `invasion_date' - 30
local invasion_ul = `invasion_date' + 30 

* Face-to-face countries 
tab cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=.

* Self-completion countries
tab cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != . 

* Capturing countries with overlap for next analyses 
levelsof cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=. & cntry != "Iceland", local(cntries_face)

levelsof cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != ., local(cntries_self) 


*** Power analysis by country based on key variable implvdm


foreach c of local cntries_face{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_face < `invasion_date' 
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'" 
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'" 
		replace power_treatment = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"  
		replace power_treatment = . if interview_face == `cut' & cntry == "`c'"  
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Droping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_face{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Switzerland.gph"	///
			"${aux}power_Greece.gph"	///
			"${aux}power_Italy.gph"	///
			"${aux}power_Montenegro.gph"	///
			"${aux}power_Macedonia.gph"	///
			"${aux}power_Netherlands.gph"	///
			"${aux}power_Norway.gph"	///
			"${aux}power_Portugal.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 				
				
* graph export "${output}power_face_to_face_alternative.eps", replace 

* Returning to master frame
frame change default 
	

*** Treatment variables based on power analyses (by country)


* Treatment variable for power based on half-standard deviation change 
gen treatment_halfsd = . 
foreach c of local cntries_face{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'"
	replace treatment_halfsd = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"
	replace treatment_halfsd = . if interview_face == `cut' & cntry == "`c'"
	if `bandwidth_half_`c'' == 0{
		replace treatment_halfsd = . if cntry == "`c'"
	}
}

* Power analysis by country based on key variable implvdm


foreach c of local cntries_self{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_selcomp < `invasion_date' 
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'" 
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'" 
		replace power_treatment = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"  
		replace power_treatment = . if interview_selcomp == `cut' & cntry == "`c'"  
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Droping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_self{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Spain.gph"	///
			"${aux}power_Poland.gph"	///
			"${aux}power_Serbia.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 
				
* graph export "${output}power_selfcompletion_alternative.eps", replace 

* Returning to master frame
frame change default 
	

* Treatment variables based on power analyses (by country)

* Treatment variable for power based on half-standard deviation change 
foreach c of local cntries_self{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'"
	replace treatment_halfsd = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"
	replace treatment_halfsd = . if interview_selcomp == `cut' & cntry == "`c'"
	if `bandwidth_half_`c''  == 0{
		replace treatment_halfsd = . if cntry == "`c'"
	}
}

*** Minor corrections 

foreach var of varlist treatment_*{
	replace `var' = . if cntry == "Norway"
}


* Bandwidth by country and variable
levelsof cntry if treatment_halfsd != ., local(countries)

foreach x of local countries{
	frame change observations_`x'
	sort Bandwidth
	sum Bandwidth if half_sd == 1
	local half_`x' = r(min)
	sum Bandwidth if third_sd == 1
	local third_`x' = r(min)
} 

frame change default

** Locals for plots labels with country name + N + Bandwidth

levelsof cntry if treatment_halfsd != ., local(countries)
foreach x of local countries{
	local band = `half_`x''
	sum treatment_halfsd if cntry == "`x'"
	local observations = r(N)
	local third_lab_`x' = `"`x' = `""`x'" "N = `observations'" "Bandwidth = `band' days""'"'
	local band_plot_`x' = `band'
}


*** Analysis

matrix drop _all

* Matrices to store results 
levelsof cntry if treatment_halfsd != ., local(countries)
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
	
	matrix `var'_corr = J(1,10,.)
	matrix colnames `var'_corr = `countries'
}

* Analysis
local i = 0
foreach x of local countries{
	
	local ++ i 
	
	* Balancing weight
	if "`x'" != "Montenegro"{
	ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							if cntry == "`x'",  generate(balance) targets(2) tolerance(.050)
	
	svyset [pweight=balance]
	}
	* Analysis 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'", level(95)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'", level(95)
		}		
		local b = r(table)["b", 1]
		local ll = r(table)["ll", 1]
		local ul = r(table)["ul", 1]
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'", level(99)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'", level(99)
		}	
		local ll99 = r(table)["ll", 1]
		local ul99 = r(table)["ul", 1]
		
	* Storing results
	matrix `var'[1, `i'] = `b' 
	matrix `var'_CI[1, `i'] = `ll' 
	matrix `var'_CI[2, `i'] = `ul' 
	matrix `var'_CI[3, `i'] = `ll99' 
	matrix `var'_CI[4, `i'] = `ul99'
	}
	
	* Multiple hypothesis testing corrections 
	if "`x'" == "Greece"{
		local num = 9
	} 
	if "`x'" == "Italy"{
		local num = 12
	} 
	if "`x'" == "Macedonia"{
		local num = 14
	} 
	if "`x'" == "Montenegro"{
		local num = 15
	} 
	if "`x'" == "Netherlands"{
		local num = 16
	} 
	if "`x'" == "Poland"{
		local num = 18
	} 
	if "`x'" == "Portugal"{
		local num = 19
	} 
	if "`x'" == "Serbia"{
		local num = 20
	} 
	if "`x'" == "Spain"{
		local num = 23
	} 
	if "`x'" == "Switzerland"{
		local num = 25
	} 	
	
	if "`x'" != "Montenegro"{
		rwolf2 	(svy: reg implvdm treatment_halfsd  if cntrynum == `num')	///
				(svy: reg accalaw treatment_halfsd  if cntrynum == `num') ///
				(svy: reg stfdem treatment_halfsd  if cntrynum == `num')	///
				(svy: reg fairelc treatment_halfsd  if cntrynum == `num') ///
				(svy: reg rghmgpr treatment_halfsd  if cntrynum == `num')	///
				(svy: reg votedir treatment_halfsd  if cntrynum == `num')	///
				(svy: reg cttresa treatment_halfsd  if cntrynum == `num')	///
				(svy: reg gvctzpv treatment_halfsd  if cntrynum == `num')	///
				(svy: reg grdfinc treatment_halfsd  if cntrynum == `num')	///
				(svy: reg viepol treatment_halfsd  if cntrynum == `num')	///
				(svy: reg wpestop treatment_halfsd  if cntrynum == `num')	///
				(svy: reg keydec treatment_halfsd  if cntrynum == `num')	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
		if "`x'" == "Montenegro"{
		rwolf2 	(reg implvdm  treatment_halfsd if cntrynum == `num')	///
				(reg accalaw  treatment_halfsd if cntrynum == `num') ///
				(reg stfdem  treatment_halfsd if cntrynum == `num')	///
				(reg fairelc  treatment_halfsd if cntrynum == `num') ///
				(reg rghmgpr  treatment_halfsd if cntrynum == `num')	///
				(reg votedir  treatment_halfsd if cntrynum == `num')	///
				(reg cttresa  treatment_halfsd if cntrynum == `num')	///
				(reg gvctzpv  treatment_halfsd if cntrynum == `num')	///
				(reg grdfinc  treatment_halfsd if cntrynum == `num')	///
				(reg viepol  treatment_halfsd  if cntrynum == `num')	///
				(reg wpestop  treatment_halfsd if cntrynum == `num')	///
				(reg keydec  treatment_halfsd  if cntrynum == `num')	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
	* Storing resuls multiple hypotheses corrections
	matrix multiple =  e(RW) 
	local r = 0 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		local ++ r 
		matrix `var'_corr[1, `i'] = multiple[`r',3] 
	}
	
	
	* Drop balancing weights variable
	if "`x'" != "Montenegro"{
			drop balance
	}
}

* Combined plot 
* Regular plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[3] `var'_CI[4]) (`var'_CI[1] `var'_CI[2])))"
}

* Corrections plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	matrix `var'_corr = `var'_corr
	local coefplot_`var'_corr = "(matrix(`var'_corr), label(`label'))"
	matrix colnames `var' = `countries'
}

* Plot for first group of variables
coefplot	`coefplot_implvdm' `coefplot_accalaw' `coefplot_stfdem' `coefplot_fairelc' 	`coefplot_rghmgpr' `coefplot_cttresa' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_implvdm_corr' `coefplot_accalaw_corr' `coefplot_stfdem_corr' `coefplot_fairelc_corr' 	`coefplot_rghmgpr_corr' `coefplot_cttresa_corr'	///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law"))	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))
				
graph export "${output}small_effects_var1_alternative.eps", replace 
			
* Plot for second group of variables
coefplot	`coefplot_votedir'  `coefplot_gvctzpv' `coefplot_grdfinc' `coefplot_viepol'	`coefplot_wpestop' `coefplot_keydec' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_votedir_corr'  `coefplot_gvctzpv_corr' `coefplot_grdfinc_corr' `coefplot_viepol_corr'	`coefplot_wpestop_corr' `coefplot_keydec_corr' ///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))

graph export "${output}small_effects_var2_alternative_date.eps"	, replace 	


*** Analysis: POOLED SAMPLE SMALL EFFECTS

* Matrices to store results 
matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balacing weight 
ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis 
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' treatment_halfsd i.cntrynum, level(95)
	local b = r(table)["b", 1]
	local ll = r(table)["ll", 1]
	local ul = r(table)["ul", 1]
	svy: reg `var' treatment_halfsd i.cntrynum, level(99)
	local ll99 = r(table)["ll", 1]
	local ul99 = r(table)["ul", 1]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}

* Analsysis multiple hypothesis test corrections
rwolf2 	(svy: reg implvdm treatment_halfsd i.cntrynum)	///
		(svy: reg accalaw treatment_halfsd i.cntrynum) ///
		(svy: reg stfdem treatment_halfsd i.cntrynum)	///
		(svy: reg fairelc treatment_halfsd i.cntrynum) ///
		(svy: reg rghmgpr treatment_halfsd i.cntrynum)	///
		(svy: reg votedir treatment_halfsd i.cntrynum)	///
		(svy: reg cttresa treatment_halfsd i.cntrynum)	///
		(svy: reg gvctzpv treatment_halfsd i.cntrynum)	///
		(svy: reg grdfinc treatment_halfsd i.cntrynum)	///
		(svy: reg viepol treatment_halfsd i.cntrynum)	///
		(svy: reg wpestop treatment_halfsd i.cntrynum)	///
		(svy: reg keydec treatment_halfsd i.cntrynum)	///
		, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
		reps(`replications') verbose seed(`seed_corrections')

matrix multiple =  e(RW) 		
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	matrix results_corr[1,`c'] = multiple[`c',3] 
}



drop balance

* Combined plot 
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			, bylabel((A) Estimation)	///
			|| (matrix(results_corr)) ///
			, bylabel((B) Multiple testing corrections)	///
			||,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	///
			byopts(xrescale)

addplot 1: ,	legend(off) 	///
				xlabel(-1(1)1) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(off) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))


graph export "${output}small_effects_pooled_alternative_date.eps", replace


*****************************************
* Excluding paper-and-pencil interviews
*****************************************
* Defining new date
matrix drop _all

levelsof cntry if treatment_halfsd != ., local(countries)
foreach x of local countries{
	frame drop observations_`x'
}

frame drop observations_Norway

drop treatment_*

* Locals for date and initial bandwidth arround invasion date (1 month)
local invasion_date = date("24feb2022","DMY")
local invasion_ll = `invasion_date' - 30
local invasion_ul = `invasion_date' + 30 

* Face-to-face countries 
tab cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=.

* Self-completion countries
tab cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != . 

* Capturing countries with overlap for next analyses 
levelsof cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=. & cntry != "Iceland", local(cntries_face)

levelsof cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != ., local(cntries_self) 


*** Power analysis by country based on key variable implvdm


foreach c of local cntries_face{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_face < `invasion_date' & mode != 4
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'" 
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'"  & mode != 4
		replace power_treatment = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'"  & mode != 4
		replace power_treatment = . if interview_face == `cut' & cntry == "`c'"  & mode != 4
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" & mode != 4 // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 & mode != 4 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 & mode != 4 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Droping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_face{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Switzerland.gph"	///
			"${aux}power_Greece.gph"	///
			"${aux}power_Italy.gph"	///
			"${aux}power_Montenegro.gph"	///
			"${aux}power_Macedonia.gph"	///
			"${aux}power_Netherlands.gph"	///
			"${aux}power_Norway.gph"	///
			"${aux}power_Portugal.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 				
				
* graph export "${output}power_face_to_face_alternative.eps", replace 

* Returning to master frame
frame change default 
	

*** Treatment variables based on power analyses (by country)


* Treatment variable for power based on half-standard deviation change 
gen treatment_halfsd = . 
foreach c of local cntries_face{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_face < `cut' & interview_face > `low' & cntry == "`c'" & mode != 4
	replace treatment_halfsd = 1 if interview_face >= `cut' & interview_face <= `up' & cntry == "`c'" & mode != 4
	replace treatment_halfsd = . if interview_face == `cut' & cntry == "`c'" & mode != 4
	if `bandwidth_half_`c'' == 0{
		replace treatment_halfsd = . if cntry == "`c'" 
	}
}

* Power analysis by country based on key variable implvdm


foreach c of local cntries_self{
	
	* Parameters for power calculation. Mean and SD of main variable of interest
	sum implvdm if cntry == "`c'" & interview_selcomp < `invasion_date' & mode != 4
	local mean = r(mean)
	local sd = r(sd) 
	local mitad = `sd' / 2 
	local tercio = `sd' / 3 
	local half = `mean' - `mitad'
	local third = `mean' - `tercio'
	
	* Definition of matrix to store the results
	matrix results_`c' = J(30,6,.)
	
	* Loop to calculate power for each bandwidth
	forval x = 1/30{
		
		* Treatment variable definition
		local days_after = `x' // Bandwidth 
		gen power_treatment = . if cntry == "`c'"  & mode != 4
		local cut = `invasion_date'
		local up = `cut' + `days_after' 
		local low = (`cut' - `days_after') - 1 
		replace power_treatment = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'" & mode != 4
		replace power_treatment = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'"  & mode != 4
		replace power_treatment = . if interview_selcomp == `cut' & cntry == "`c'"  & mode != 4
		
			
		* Geting the number of observations
		sum power_treatment if cntry == "`c'" & mode != 4 // Total 
		local total = r(N)		
		sum power_treatment if cntry == "`c'" & power_treatment == 0 & mode != 4 // Control 
		local total_control = r(N)
		sum power_treatment if cntry == "`c'" & power_treatment == 1 & mode != 4 // Treatment 
		local total_treatment = r(N)
		
		* Estimating power for each bandwidth
		capture noisily power twomeans `mean' `third', sd(`sd') n1(`total_control') n2(`total_treatment')  
		local pow_third = r(power)
		capture noisily power twomeans `mean' `half', sd(`sd') n1(`total_control') n2(`total_treatment') 
		local pow_half = r(power)
		
		
		* Inputing matrix
		matrix results_`c'[`x', 1] = `x',`total',`total_control',`total_treatment',`pow_third',`pow_half'
		matrix colnames results_`c' =  "Bandwidth" "Observations" "Obs_control" "Obs_treatment" "Power_third_SD" "Power_half_SD"
		
		* Droping variables 
		drop power_treatment
	}
}

* Creating plots and variable capturing the optimal bandwidth by country
foreach c of local cntries_self{
	
	* Creating frame to store and plot the results
	frame create observations_`c'
	frame change observations_`c'
	
	* Transfering matrix to frame
	svmat results_`c', names(col)
	
	* Calculating parameters for plot
	sum Observations
	local max = r(max)
	local steps = r(max)/7
	
	* Generating plot
	graph twoway	(bar Observations Bandwidth, fcolor(gs4) ///
					yaxis(2) ylabel(0(`steps')`max', format(%6.1g) labsize(vsmall) axis(2)) ytitle("") ///
					yscale(axis(2) alt)) ///
					(rbar Observations Obs_treatment Bandwidth, fcolor(gs13) lcolor(gs1)  lwidth(vthin) yaxis(2)) ///
					(function y=.8, ra(0 31) lpattern(solid) lcolor(red) lwidth(vthin)) ///
					(line Power_half_SD Bandwidth, lcolor(blue) lpattern(longdash) yaxis(1) yscale(axis(1) alt)) ///
					(line Power_third_SD Bandwidth, lcolor(orange) lpattern(shortdash)  yaxis(1) ///
					ytitle("") ylabel(0(.2)1, labsize(vsmall) gmin gmax axis(1)) ///
					xtitle("Bandiwdth (+/- days) around invasion") xlabel(0(5)30, labsize(vsmall)) ///
					legend(ti("Effect size", size(small)) order(4 "1/2 SD change" 5 "1/3 SD change") size(small) rows(1) pos(6))	///
					ti("`c'", size(small)))
	
	* Saving plot
	graph save "${aux}power_`c'.gph", replace
	
	* Generating variables 
	gen half_sd = 1 if Power_half_SD >= 0.8
	gen third_sd = 1 if Power_third_SD >= 0.8
	sum Bandwidth if half_sd == 1
	if `r(N)' > 0{
		local bandwidth_half_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_half_`c' = 0
	}
	sum Bandwidth if third_sd == 1
	if `r(N)' > 0{
		local bandwidth_third_`c' = r(min)
	}
	if `r(N)' == 0{
		local bandwidth_third_`c' = 0
	}
}

grc1leg2	"${aux}power_Spain.gph"	///
			"${aux}power_Poland.gph"	///
			"${aux}power_Serbia.gph"	///
			,	l1("Statistical power", size(small)) ///
				r1("Number of cases in treatment and control groups", size(small)) 
				
* graph export "${output}power_selfcompletion_alternative.eps", replace 

* Returning to master frame
frame change default 
	

* Treatment variables based on power analyses (by country)

* Treatment variable for power based on half-standard deviation change 
foreach c of local cntries_self{
	local cut = `invasion_date'
	local up = `cut' + `bandwidth_half_`c'' 
	local low = (`cut' - `bandwidth_half_`c'') - 1 
	replace treatment_halfsd = 0 if interview_selcomp < `cut' & interview_selcomp > `low' & cntry == "`c'" & mode != 4
	replace treatment_halfsd = 1 if interview_selcomp >= `cut' & interview_selcomp <= `up' & cntry == "`c'" & mode != 4
	replace treatment_halfsd = . if interview_selcomp == `cut' & cntry == "`c'" & mode != 4
	if `bandwidth_half_`c''  == 0{
		replace treatment_halfsd = . if cntry == "`c'" & mode != 4
	}
}

*** Minor corrections 

foreach var of varlist treatment_*{
	replace `var' = . if cntry == "Norway"
}


* Bandwidth by country and variable
levelsof cntry if treatment_halfsd != ., local(countries)

foreach x of local countries{
	frame change observations_`x'
	sort Bandwidth
	sum Bandwidth if half_sd == 1
	local half_`x' = r(min)
	sum Bandwidth if third_sd == 1
	local third_`x' = r(min)
} 

frame change default

** Locals for plots labels with country name + N + Bandwidth

levelsof cntry if treatment_halfsd != ., local(countries)
foreach x of local countries{
	local band = `half_`x''
	sum treatment_halfsd if cntry == "`x'"
	local observations = r(N)
	local third_lab_`x' = `"`x' = `""`x'" "N = `observations'" "Bandwidth = `band' days""'"'
	local band_plot_`x' = `band'
}


*** Analysis

matrix drop _all

* Matrices to store results 
levelsof cntry if treatment_halfsd != ., local(countries)
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	matrix `var' = J(1,10,.)
	matrix `var'_CI = J(4,10,.)
	matrix colnames `var' = `countries'
	matrix colnames `var'_CI = `countries'
	
	matrix `var'_corr = J(1,10,.)
	matrix colnames `var'_corr = `countries'
}

* Analysis
local i = 0
foreach x of local countries{
	
	local ++ i 
	
	* Balancing weight
	if "`x'" != "Montenegro"{
	ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
							if cntry == "`x'" & mode != 4,  generate(balance) targets(2) tolerance(.050)
	
	svyset [pweight=balance]
	}
	* Analysis 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'" & mode != 4, level(95)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'" & mode != 4, level(95)
		}		
		local b = r(table)["b", 1]
		local ll = r(table)["ll", 1]
		local ul = r(table)["ul", 1]
		if "`x'" != "Montenegro"{
			quietly: svy: reg `var' treatment_halfsd  if cntry == "`x'" & mode != 4, level(99)
		}
		if "`x'" == "Montenegro"{
			quietly: reg `var' treatment_halfsd if cntry == "`x'" & mode != 4, level(99)
		}	
		local ll99 = r(table)["ll", 1]
		local ul99 = r(table)["ul", 1]
		
	* Storing results
	matrix `var'[1, `i'] = `b' 
	matrix `var'_CI[1, `i'] = `ll' 
	matrix `var'_CI[2, `i'] = `ul' 
	matrix `var'_CI[3, `i'] = `ll99' 
	matrix `var'_CI[4, `i'] = `ul99'
	}
	
	* Multiple hypothesis testing corrections 
	if "`x'" == "Greece"{
		local num = 9
	} 
	if "`x'" == "Italy"{
		local num = 12
	} 
	if "`x'" == "Macedonia"{
		local num = 14
	} 
	if "`x'" == "Montenegro"{
		local num = 15
	} 
	if "`x'" == "Netherlands"{
		local num = 16
	} 
	if "`x'" == "Poland"{
		local num = 18
	} 
	if "`x'" == "Portugal"{
		local num = 19
	} 
	if "`x'" == "Serbia"{
		local num = 20
	} 
	if "`x'" == "Spain"{
		local num = 23
	} 
	if "`x'" == "Switzerland"{
		local num = 25
	} 	
	
	if "`x'" != "Montenegro"{
		rwolf2 	(svy: reg implvdm treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg accalaw treatment_halfsd  if cntrynum == `num' & mode != 4) ///
				(svy: reg stfdem treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg fairelc treatment_halfsd  if cntrynum == `num' & mode != 4) ///
				(svy: reg rghmgpr treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg votedir treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg cttresa treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg gvctzpv treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg grdfinc treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg viepol treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg wpestop treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(svy: reg keydec treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
		if "`x'" == "Montenegro"{
		rwolf2 	(reg implvdm  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg accalaw  treatment_halfsd if cntrynum == `num' & mode != 4) ///
				(reg stfdem  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg fairelc  treatment_halfsd if cntrynum == `num' & mode != 4) ///
				(reg rghmgpr  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg votedir  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg cttresa  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg gvctzpv  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg grdfinc  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg viepol  treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				(reg wpestop  treatment_halfsd if cntrynum == `num' & mode != 4)	///
				(reg keydec  treatment_halfsd  if cntrynum == `num' & mode != 4)	///
				, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
				reps(`replications') verbose seed(`seed_corrections')
	}
	
	* Storing resuls multiple hypotheses corrections
	matrix multiple =  e(RW) 
	local r = 0 
	foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
		local ++ r 
		matrix `var'_corr[1, `i'] = multiple[`r',3] 
	}
	
	
	* Drop balancing weights variable
	if "`x'" != "Montenegro"{
			drop balance
	}
}

* Combined plot 
* Regular plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	local coefplot_`var' = "(matrix(`var'), label(`label')  ci((`var'_CI[3] `var'_CI[4]) (`var'_CI[1] `var'_CI[2])))"
}

* Corrections plots specification 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local label : variable label `var' 
	matrix `var'_corr = `var'_corr
	local coefplot_`var'_corr = "(matrix(`var'_corr), label(`label'))"
	matrix colnames `var' = `countries'
}

* Plot for first group of variables
coefplot	`coefplot_implvdm' `coefplot_accalaw' `coefplot_stfdem' `coefplot_fairelc' 	`coefplot_rghmgpr' `coefplot_cttresa' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_implvdm_corr' `coefplot_accalaw_corr' `coefplot_stfdem_corr' `coefplot_fairelc_corr' 	`coefplot_rghmgpr_corr' `coefplot_cttresa_corr'	///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Importance democracy" 6 "Strong leader above law" 9 "Free elections" 12 "Satisfaction democracy" 15 "Minority rights" 18 "Rule of law"))	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))
				
graph export "${output}paper_pencli_excluded_var1.eps", replace 
			
* Plot for second group of variables
coefplot	`coefplot_votedir'  `coefplot_gvctzpv' `coefplot_grdfinc' `coefplot_viepol'	`coefplot_wpestop' `coefplot_keydec' ///
			, bylabel((A) Estimation)	///
			|| `coefplot_votedir_corr'  `coefplot_gvctzpv_corr' `coefplot_grdfinc_corr' `coefplot_viepol_corr'	`coefplot_wpestop_corr' `coefplot_keydec_corr' ///
			, bylabel((B) Multiple testing corrections)	///
			|| ,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall))	///
			coeflabels(`third_lab_Greece' `third_lab_Italy' `third_lab_Macedonia' `third_lab_Montenegro' `third_lab_Netherlands' `third_lab_Poland' `third_lab_Portugal' `third_lab_Serbia' `third_lab_Spain' `third_lab_Switzerland')	///
			 legend(cols(3) rows(2)) byopts(xrescale)

addplot 1: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				xlabel(-3(1)3) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(order(3 "Referenda" 6 "Protection poverty" 9 "Red. income differences" 12 "Views people prevail over elite" 15 "Will of people cannot be stopped" 18 "Key decisions national not EU")) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))

graph export "${output}paper_pencli_excluded_var2.eps"	, replace 	


*** Analysis: POOLED SAMPLE SMALL EFFECTS

* Matrices to store results 
matrix drop _all

matrix results = J(1,12,.)
matrix results_CI = J(4,12,.)
matrix colnames results = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec
matrix colnames results_CI = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

matrix results_corr = J(1,12,.)
matrix colnames results_corr = implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec

* Balacing weight 
ebalance treatment_halfsd	eduyrs gndr agea active	voted						///
						if mode !=4	,  generate(balance) targets(3) tolerance(.050)
svyset [pweight=balance]

* Analysis 
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	* Estimation
	svy: reg `var' treatment_halfsd i.cntrynum if mode !=4	, level(95)
	local b = r(table)["b", 1]
	local ll = r(table)["ll", 1]
	local ul = r(table)["ul", 1]
	svy: reg `var' treatment_halfsd i.cntrynum if mode !=4	, level(99)
	local ll99 = r(table)["ll", 1]
	local ul99 = r(table)["ul", 1]	
	
	* Storing results
	matrix results[1,`c'] = `b' 
	matrix results_CI[1,`c'] = `ll' 
	matrix results_CI[2,`c'] = `ul' 
	matrix results_CI[3,`c'] = `ll99' 
	matrix results_CI[4,`c'] = `ul99'
}

* Analsysis multiple hypothesis test corrections
rwolf2 	(svy: reg implvdm treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg accalaw treatment_halfsd i.cntrynum if mode !=4	) ///
		(svy: reg stfdem treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg fairelc treatment_halfsd i.cntrynum if mode !=4	) ///
		(svy: reg rghmgpr treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg votedir treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg cttresa treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg gvctzpv treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg grdfinc treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg viepol treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg wpestop treatment_halfsd i.cntrynum if mode !=4	)	///
		(svy: reg keydec treatment_halfsd i.cntrynum if mode !=4	)	///
		, indepvars(treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd, treatment_halfsd)	///
		reps(`replications') verbose seed(`seed_corrections')

matrix multiple =  e(RW) 		
local c 0 
foreach var of varlist implvdm accalaw stfdem fairelc  rghmgpr votedir cttresa gvctzpv grdfinc viepol wpestop keydec{
	local ++ c
	matrix results_corr[1,`c'] = multiple[`c',3] 
}

drop balance

* Combined plot 
coefplot	(matrix(results), ci((results_CI[3] results_CI[4]) (results_CI[1] results_CI[2]))) ///
			, bylabel((A) Estimation)	///
			|| (matrix(results_corr)) ///
			, bylabel((B) Multiple testing corrections)	///
			||,  msize(vsmall)  ciopts(lwidth(thin medium)) ///
			ylabel(, labsize(vsmall) notick) legend(size(vsmall)) legend(off) 	///
			mlabel format(%9.2g) mlabposition(12) mlabgap(*1.5) mlabsize(vsmall)	///
			byopts(xrescale)

addplot 1: ,	legend(off) 	///
				xlabel(-1(1)1) xline(0, lpattern(solid)) norescaling ///
				b1title("Point estimate and confidence interval", size(vsmall))

addplot 2: ,	legend(off) 	///
				norescaling ///
				xlabel(0(.1)1) xline(0.05, lpattern(dash) lcolor(red)) xline(0.01, lpattern(dash) lcolor(green))  	/// 
				b1title("Romano-Wolf stepdown p-values for multiple hypothesis testing", size(vsmall))


graph export "${output}paper_and_pencil_excluded_pooled.eps", replace

************************************
* Benchmarking of statistical power 
************************************

* Locals for date and initial bandwidth arround invasion date (1 month)
local invasion_date = date("24feb2022","DMY")
local invasion_ll = `invasion_date' - 30
local invasion_ul = `invasion_date' + 30 

* Face-to-face countries 
tab cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=.

* Self-completion countries
tab cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != . 

* Capturing countries with overlap for next analyses 
levelsof cntry if interview_face >= `invasion_ll' & interview_face <= `invasion_ul' & interview_face !=. & cntry != "Iceland", local(cntries_face)

levelsof cntry if interview_selcomp >= `invasion_ll' & interview_selcomp <= `invasion_ul' & interview_selcomp != ., local(cntries_self) 

* Distribution of key variable in the pooled sample 

sum implvdm if included == 1 
di "Half a standard deviation is equal to " `r(sd)'/2
di "A third of a standard deviation is equal to " `r(sd)'/3

* Comparison to other key variables of interest


** Education

* Recoding education 
gen education_level = . 
replace education_level = 1 if edulvlb <=113
replace education_level = 1 if edulvlb > 113 & edulvlb <= 229 
replace education_level = 3 if edulvlb > 229 &  edulvlb <= 323 
replace education_level = 4 if edulvlb > 323 & edulvlb <= 423
replace education_level = 5 if edulvlb > 423 & edulvlb <= 520
replace education_level = 6 if edulvlb > 520 & edulvlb <= 800
replace education_level = . if edulvlb == 5555

label def edu 1"Lower secondary (or lesss)" 3"Upper secondary" 4"Post-secondary non-tertiary" 5"Tertiary (no university)" 6"University" 

label val education_level edu

reg implvdm i.education_level if included == 1

** Political interest
reg implvdm i.polintr if included == 1


**** Generating summary table with pooled sampel and by country 

gen interview = interview_face
replace interview = interview_selcomp if interview == .

local cut = date("24feb2022","DMY")

*** Pooled 
sum implvdm if included == 1 & interview < `cut'
local half = round(`r(sd)' / 2, .01)
local third = round(`r(sd)' / 3, .01)

reg implvdm i.education_level if included == 1 & interview < `cut'
local edu = round(r(table)[1,5], .01)

reg implvdm i.polintr if included == 1 & interview < `cut'
local interest = round(r(table)[1,4], .01)


matrix pooled = `half', `third', `interest', `edu'

matrix rownames pooled = "Pooled sample"

*** By country 
levelsof cntry if included == 1, local(levels)
foreach x of local levels{
	sum implvdm if cntry == "`x'" & interview < `cut'
	local half = round(`r(sd)' / 2, .01)
	local third = round(`r(sd)' / 3, .01)

	reg implvdm i.education_level if cntry == "`x'" & interview < `cut'
	local edu = round(r(table)[1,"6.education_level"], .01)

	reg implvdm i.polintr if cntry == "`x'" & interview < `cut'
	local interest = round(r(table)[1,4], .01)

	matrix `x' = `half', `third', `interest', `edu'

	matrix rownames `x' = "`x'"
	
}

* Generating table

matrix power_benchmark = pooled \ Switzerland \ Spain \ Serbia \ Portugal \ Poland \ Netherlands \ Montenegro \ Macedonia \ Italy \ Greece

matrix colnames power_benchmark = "1/2 SD" "1/3 SD" "Political interest" "Education"

esttab matrix(power_benchmark) using "${output}power_benchmark.tex", tex replace nomtitles


* Country differences
levelsof cntry if included == 1, local(levels)
foreach x of local levels{
	di "`x'"
	sum implvdm if cntry == "`x'"
}
		



